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ABSTRACT 

A method for integrating the chemical equations associated with nuclear com- 
bustion at high temperature is presented and extensively checked. Following the 
idea of Miiller (1986), the feedback between nuclear rates and temperature was 
taken into account by simultaneously computing molar fraction changes and tem- 
perature response in the same matrix. The resulting algorithm is very stable and 
efficient at calculating nuclear combustion in explosive scenarios, especially in 
those situations where the reacting material manages to climb to the nuclear 
statistical equilibrium regime. The numerical scheme may be useful not only for 
those who carry out hydrodynamical simulations of explosive events, but also as 
a tool to investigate the properties of a nuclear system approaching equilibrium 
through a variety of thermodynamical trajectories. 

Subject headings: nuclear reaction networks, explosive nucleosynthesis- super- 
novae: general 



1. Introduction 

Explosive phenomena in stars are often associated with the rapid release of nuclear en- 
ergy and a lavish production of nuclear species. In any calculation dealing with novae or 
supernovae explosions it is highly desirable to implement the nuclear part of hydrodynamic 
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codes as better as possible in order to get a suitable representation of these explosions and 
nucleosynthetic yields. A complication arises because of the huge quantity of nuclei involved, 
which makes it difficult to simultaneously handle hydrodynamics and nuclear combustion. 
Recently, however, there have been performed hydrodynamic simulations in spherical sym- 
metry incorporating several thousand nuclei (Rauscher et al. 2001). Also, calculations in 
two or three dimensions of both kinds of supernovae are quite common. In these multidi- 
mensional calculations the mission of the nuclear network is to provide a reasonable nuclear 
energy rate at a minimum computational cost. This is usually accomplished by using small 
networks of 7 or 13 isotopes (Timmes, Hoffman & Woosley 2000). Detailed nucleosynthesis 
is calculatedted a posteriori by using a postprocessing technique with a larger network. 

The most common procedure used for integrating the nuclear networks handles the 
discretized chemical equations in their implicit form. This enhances the stability of the 
nuclear system, which is very rigid (Arnett 1996), allowing it to take longer time-steps. It is 
worthy of note that despite the great dependence of the rates on the temperature, the effect 
of the latter is rarely incorporated implicitly into numerical codes. This enforces to make a 
reliable extrapolation of the temperature and to take small time-steps to guarantee stability, 
especially when the nuclei approach the nuclear statistical equilibrium (NSE) regime. Several 
methods have been suggested for handling nuclear combustion around the NSE regime. The 
most simple and fastest way is to use a conventional integration method when the nuclear 
system is not too close to equilibrium and to rely on statistical techniques to treat the phase 
in which direct and reverse reactions are almost or completely balanced. Nevertheless, this 
leads to the problem of how to make a smooth transition between the two regimes. Another 
possibility is to always use a standard integration method, which takes care of adequately 
restricting the time-step near equilibrium. In this case the computational effort increases. 
There are also mixed schemes (Hix & Thielemann 1996) in which reduced networks are used 
to describe the evolution of light and a few intermediate-mass elements (generally below 
28 Si or 32 S) and a statistical approach is taken for handling the remaining nuclei. While 
this method is a very promising tool for explosive nucleosynthesis studies, due to its ability 
to handle hundreds of nuclei with a low computational burden, the decoupling between the 
thermal evolution and nuclear kinetics again imposes a constraint on time-steps near NSE. 

The importance of solving the nuclear network by coupling molar fractions and tem- 
perature in the same matrix was demonstrated by Miiller(1986), who used a small network 
of 13 nuclei, from helium to nickel. He showed how the inclusion of the temperature effects 
was decisive in avoiding the numerical oscillations which otherwise appear during and after 
relaxation to the NSE regime. However, in that paper there were presented only very few test 
calculations, limited to isochoric combustion. Moreover, the use of an entropy equation in- 
stead of the most frequently used energy equation was preferred to calculate the temperature 
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evolution. 

In this paper we extend Miiller's previous work to consider a larger nuclear network and 
to analyze in greater detail the transition between time-dependent combustion to complete 
nuclear equilibrium and viceversa. We also propose a method for integrating based on the 
energy equation that is robust and straightforward to implement, giving very good results in 
all the realized tests. In addition we have exhaustively checked our numerical scheme against 
a large variety of thermodynamic conditions of astrophysical interest. These include not only 
isochoric combustion, as considered in Miiller seminal work, bul also adiabatic, isothermal 
and isobaric processes, as well as the particular case of combustion within a strongly coupled 
coulombic plasma. 

Nowadays the effort to write the proposed numerical scheme for larger networks (even of 
'kilo-nuclide' size) is considerably less than a decade ago, as we have the great advantage of 
a recently published library of theoretical reactions rates (i.e. Rauscher & Thielemann 2000 
or the NACRE compilation, Angulo et al. 1999). An efficient algorithm of this kind could 
be suitable for handling a variety of vigorous nuclear combustion modes such as detonation, 
flames or silicon photodisintegration. It is especially useful whenever the operator splitting 
technique has to be applied to save computational time. Currently, this often happens 
in hydrodynamic calculations of both kinds of supernovae; in these situations the implicit 
thermal coupling could help to enlarge the time-steps without degrading neither the accuracy 
nor the stability of the system. There is also the additional advantage to this scheme that it 
is more compatible with the philosophy behind the post-processing technique. A necessary 
condition for obtaining a detailed nucleosynthesis from the thermal history of an explosive 
event, calculated beforehand using a reduced network, is that the small and large networks 
give a similar energy release. Integrating the chemical equations along with the energy 
conservation helps to fulfill this requirement. 

In § 2 we describe the method of calculation and its implementation in two nuclear 
networks with 14 nuclei and 86 nuclei respectively. The ability of the scheme to monitor the 
quasi (QSE) and complete statistical equilibria is discussed. In § 3 the results of three tests 
are provided: adiabatic expansion from the NSE stage, hydrostatic and explosive silicon 
burning and nuclear flame propagation, all of them of evident interest to current studies 
on supernovae. The role of Coulomb corrections on the chemical composition of burned 
matter at high densities is also discussed in § 3. Finally, § 4 is devoted to summarizing the 
conclusions. 
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2. Method used for integrating the nuclear network 



To integrate the set of differential equations associated with the nuclear species we have 
chosen a scheme based on the simple and well checked implicit method by Arnett & Truran 
(1969). The main motivation behind such election is that Arnett&Truran method has been 
traditionally used in nucleosynthetic calculations, always giving satisfactory results with low 
computational effort. In addition it is straightforward to modify in order to incorporate the 
thermal feedback and other proposed changes into the basic scheme. There are some dif- 
ferences between our integration method and that used in Muller(1986) worth to mention. 
In Midler's work the evolution of the system was followed by coupling the nuclear evolu- 
tion with the entropy changes at constant density. Moreover, it was assumed that entropy 
changes were only due to the release of nuclear energy, thus no other sources of entropy were 
considered. Instead of the entropy equation our formulation uses energy conservation which 
is more often found in current calculations of explosive combustion. Moreover, the proposed 
algorithm is quite general, not restricted to isochoric processes, and potential entropy vari- 
ations coming from heat exchange were also included. Minor effects, such as the incidence 
of coulombic corrections to the ionic component of the plasma and screening enhancement 
factors have been also incorporated to our scheme. A simplified description of our method 
can be found in Garcfa-Senz&Cabezon (2003). 

Starting from a typical differential equation governing the evolution of the molar frac- 
tions of the i-specie, Yf 



where = pNx < cr,v >ij stands for particle reactions and A« is for photodisintegrations 
and decays. As usual, reaction rates between identical particles have to be divided by the 
corresponding factor (two or six for binary and ternary reactions respectively). The triple- 
alpha reaction was also incorporated to equations (1). Other possible three body reactions 
can be included in the same way. This set of equations is then linearized by taking 



where Y^ is the molar fraction of the i-nuclei, n the integration step and < 6 < 1 a 
parameter which allows the quality of the integration to be controlled. Choosing 6 = 1 leads 
to a totally implicit first order integration; (9 = 1/2 is centered (second order) implicit and 
= is totally explicit. Then, neglecting the terms in AY{ ■ AY}, the linearized version 




(1) 



yn+0 = yn + 0^ 



(2) 
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of the chemical equations is written. With the same level of accuracy we propose to take 
rpn+e _ rpn _|_ Q^rp a pp rox i ma t e the nuclear reaction rates by taking 

r,r^r,(r) + -^ OAT (3) 

If we also neglect second order terms in temperature, such as Al^AT, we can approximate 
a typical element on the right part of chemical equations (1) by: 



yn+l yn+1 r ..(jm+l) „ [yn r ..(T n )\ 9AYj + [Y? r^(T")] 9AY t + 



yn yn U 

* j dT 



9 AT + Y l n Y n r ij (T n ) (4) 



then, equation (1) reads: 



AYj 
At 



+ r k iYrOAY k + niY^AYt - I E r i3 Yf j 9AY % - E r^A^ + 

k,l k,l \ j / j 

E KJAY m - XiOAYi + \Y, r' kl Y k n Yr - E r^Y? + E - X'^A 9 AT + 

m \ k,l j m / 

E r <^ ny " - E r ^ ni ? + E ^ - w ( 5 ) 

Note that the derivative of nuclear rates with respect temperature, r'^ and A-, do not pose 
a practical problem since the recently published fitting formulae of Rauscher & Thielemann 
(2000) can be used to calculate them in an efficient way. Using that compilation through 
all stages of the combustion also leads to a complete internal coherence between direct and 
inverse rates when the equilibrium is approached, which is an essential condition for the 
proper handling of NSE. 

For a network of N species we have N linear equations with N + 1 unknowns. An 
additional equation, energy equation, is needed in order to close the system. For a process 
not involving particle difussion it writes: 



(6) 
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where P is pressure, U is the energy per gram of material, BE« represents the nuclear binding 
energy of a mol of the i-species and dQ can also include sources or sinks of energy other than 
nuclear. 

It is useful to write that equation by using the thermodynamical identity: 



fdP\ , 2 fdU 



thus, using; 

T n+e = T n + q/\t, the discretized form of the energy equation is easily written. 
For example, for an adiabatic process equation (6) reads: 



for an ideal gas of ions dU/dYi = (3/2)NAkT, which is usually much lower than the Q-value 
of the reactions. Coulomb corrections to the ideal gas must be also included when the ionic 
coupling constant parameter, Tj = 2.27x 10 5 < Z 5 / 3 > (p Y e ) 1 ^ T _1 , becomes non negligible. 
This often happens in high temperature combustion of degenerate matter because of the 
strong dependence of r, on the mean charge of the system . Such dependence also favours 
dU/dYi of the products of charged particle reactions against that of the reactant particles. 
Thus, on the whole, both contributions to dU/dYi (ideal gas and Coulomb corrections) work 
in the same direction, slightly increasing the equilibrium temperature for a given density. 
Therefore, regardless the density regime we have always kept both contributions to dU/dYi in 
equation (8). 

Interestingly the implicit treatment of temperature also allows for the proper handling of 
quasi- hydrostatic and isobaric processes. Let's suppose that the value of the actual pressure 
of the system, P, remains always close to the external pressure, P ext (t), so that P = P ext + dP. 
Then, writing dP = (dP/dT)dT + (dP/dp)dp+J2 i (dP/dY i )dY i and taken out dp the energy 
equation (eq.[6]) reads: 



-7- 



Once this equation is discretized (in a similar way that we did for eq. [6]) and solved 
jointly with equations (5) the solution, AT and AYi, can be used along with P ext to find 
Ap. Equation (9) has been used in §3.3 to solve the structure of a nuclear flame. 

From here on we will focus on the case of 6 = 1 which is expected to give the most 
stable system behaviour of the nuclear system. A recalculation of the first test presented in 
the next section with 9 = 0.5 did not lead to relevant quantitative changes in the results. 
The standard criterion to assign the integration time-step was to demand that the fractional 
change in temperature and molar fractions of the more abundant species would not exceed 
2% per model. Nevertheless, one can safely relax that criterion even to 10-15% variation in 
molar fractions without seriously modify the quantitative results. 

This set of linear equations can be efficiently solved by using standard sparse matrix 
algorithms. We in particular used the one found in Prantzos, Arnould & Arcoragi (1987). 
The resulting algorithm led to a satisfactory mass conservation, the error in the sum of mass 
fractions was always lower than 10~ 10 , close machine precision, for the calculations described 
in §3. There is not a large loss of computational efficiency when the linearized thermal 
coupling is turned on. The main source of degradation of that type is the computation for 
the derivatives of the rates rather than the increase in the size of the matrix from {N x N} to 
{(N+l) x (N+l)} . From the tests carried out we estimate that our numerical treatment leads 
to a time overload under 15% and 2% for the 14 and 86-nucleus systems respectively. Several 
tests were run which demonstrated the ability of the scheme to monitor rapid isochoric 
combustion followed by adiabatic expansion in the hydrodynamic time scale. Adding either 
diffusive terms or those related to the artificial viscosity formulation to the right side of 
equations (6) and (9) allows us to include heat transport and shock wave heating in the 
scheme. An extreme case of heat transport, the propagation of a nuclear flame, is also 
presented and checked in § 3. Even though all the test presented below were calculated 
solving the linearized equations (5) and (8) an improvement of the algorithm based on the 
iterative refinement of the solution is also given and discussed in §3.1.2. 



3. Testing the algorithm 

Three test cases were calculated in order to check the method. All of them are of an 
indisputable interest in explosive combustion: a) isochoric combustion followed by adiabatic 
expansion once complete NSE has been achieved, b) endoenergetic disintegration of silicon 
and c) isobaric nuclear flame propagation. 

The nuclear network consists of 86 nuclei spanning from 12 C to 60 Zn including a, p and 
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n particles, as depicted in Figure 1. Even though this is not a large network, but rather a 
moderate one, it is certainly more complex than the a-chain considered in Miiller (1986). 
On another note, to avoid an excessive computational burden, current multidimensional 
simulations of supernovae can only support very reduced networks. For this reason, we also 
carried out several runs of the test using an a-network of 14 nuclei, a subset of the previous 
sample of 86 nuclei. This reduced set was similar to the sets of 7 and 13 nuclei considered 
in the comparative study about the performance of small nuclear networks carried out by 
Timmes et al. (2000) . Nuclear rates on light elements ,up to 20 Ne, were taken from Caughlan 
& Fowler (1988), and from Rauscher & Thielemann (2000) above neon. 

A common feature of all the tests shown here is that nuclear burning proceeds through 
the so-called e-process, thus it manages to achieve complete nuclear equilibrium at high 
temperature. In nature, that equilibrium does not last very long because it is cut by the 
reduction in density as the expansion progresses. At some point the temperature becomes too 
low to allow the nuclear reactions to proceed and the chemical abundances freeze out. The 
precise value of the resulting abundances is a function of the initial entropy of the material 
in the NSE phase, the degree of neutronization of the material and the rate of the expansion, 
as many studies have demonstrated (see for instance Meyer, Krishnan & Clayton 1998 for 
a recent work and references therein). 

The equation of state (EOS) used in all the tests was quite complete, consisting of 
a mix of electrons treated as a partially degenerate relativistic gas, an ideal gas of ions 
with Coulomb and other minor corrections, and radiation. Although important at high 
densities, electron capture on protons and nuclei have not been included in these calculations. 
Nevertheless, they can be incorporated into the scheme with the same formulation expressed 
in equations (2) and (3). 

The treatment of electrostatic corrections deserves further discussion. At high densities 
the energy of the plasma deviates from that of a pure ideal gas of nuclei and electrons. 
Electrostatic interactions must be taken into account in this regime, even in NSE, owing 
to the strong dependence of such interactions with the average charge of the system. The 
electrostatic energy of the k-specie associated to the Coulomb interactions, Uj? (erg. g _1 ), 
can be approximated by the following expressions (Ogata & Ichimaru 1987, Yakolev & 
Shalybkov 1987): 



where Tk is the ionic coupling constant parameter of the k-nuclei, a=-0.8980, b=0. 96786, 




(10) 
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c=0.2207, d=-0.86097, (3 = 0.29561 and 7 = 1.9885. The chemical potentials, ^, can be 
easily derived from expression (10). Nevertheless, the inclusion of Coulombic interactions 
affects not only the thermodynamics, they also slightly alters the rates of electron captures 
and the rate of nuclear reactions through the corresponding screening factors. As we have 
not included weak interactions in the network we will not address the first point. The second 
point, the incorporation of screening corrections to the nuclear rates in the scheme should 
be straightforward. However, there is also necessary to compute the derivatives of such 
corrections with respect to temperature, as we did with the nuclear reaction rates. 

In the following, we are mainly interested in the behaviour of the system when it ap- 
proaches the NSE. In this case it is enough to make a simplified treatment of the screening. In 
particular, we take the enhancement factors given by Mochkovitch & Nomoto (1986) who, 
among other things, studied the equilibrium of direct and reverse reactions when strong 
screening corrections are taken into account. They found that the dominant term in the 
enhancement factor, EF, to direct capture reactions is: 



c- 



EF 



cxp 



(11) 



where A/i c is the chemical potential of the reactants minus that of the products. The above 
expression is especially adequate for non-resonant reactions at high temperature as well as 
for resonant reactions with large enough resonance energy. Note that, in this approximation 
the use of the enhancement factors given by equation (11) does not affect to the photodis- 
integration reactions calculated through the detailed balance, because the Q-value of the 
reaction becomes Q — > Qo + A/z c which exactly compensates the screening correction. In 
the limit of complete nuclear equilibrium this approach leads to the same results as the 
statistical description given by the Saha equation with the chemical potentials corrected for 
electrostactic interactions (Bravo& Garcia-Senz 1999). 



3.1. Isochoric combustion and further expansion 

Isochoric combustion is a very interesting test because the original fuel (usually helium 
or carbon and oxygen in many explosive events) goes through several modes of burning: 
initial fast burning followed by quasi-statistical and complete nuclear equilibrium. Further 
expansion will take the combustion out of NSE, which leads to a rapid freezing of the 
abundances. In hydrodynamic calculations, that regime might be associated to a detonation 
wave passing through degenerate material. 
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To begin, the combustion of a system composed of equal parts 12 C and 16 at p = 
10 9 g.cm" 3 was followed without including the implicit coupling between abundances and 
temperature. The result is shown in Figure 2, where we can see that mass fractions begin 
to oscillate as soon as the combustion approaches the equilibrium regime. This behavior is 
characteristic of any stiff system of differential equations solved with explicit rather than 
implicit schemes, demanding very short time-steps to control the system (< 10~ 9 s). The 
evolution with thermal coupling is depicted in Figure 3; here the oscillations did not show 
up and the QSE and NSE regimes were handled without problems, despite the large time- 
steps (~ 0.1 s, even larger than the hydrodynamic timescale thd, see below) taken once 
equilibrium was achieved. From Tables 1 and 2 we can see that the most abundant nuclei 
correspond to Fe-peak elements, especially 54 Fe (27%, in mass) and also to helium (almost 
20%). Such abundances and the high temperature achieved during the isochoric combustion 
are in agreement with the current picture of NSE (Arnett 1996). 

The integration scheme also allowed us to follow a rapid adiabatic expansion in the 
hydrodynamic timescale thd = 446/ y/p s immediately after the NSE regime was reached: 

p(t) = p exp[-t/T HD ] (12) 

During the expansion there was no need to impose any particular freeze-out temperature 
for the nucleonic system. In Table 1 we provide the mass fractions at five different fiducial 
temperatures achieved along the cooling path. From these numbers we can see that there 
is a large spread in the freezing temperatures. Many species are still changing even below 
T = 2 x 10 9 K, although their abundances are generally too low to be significative. The 
most abundant element, 56 Ni, freezed around T ~ 3 x 10 9 K which is not a high enough 
temperature to keep the NSE going. Thus, whenever the Saha equation is used to describe 
nuclear statistical equilibrium above 4 — 5 x 10 9 K one has to be cautious and switch to a 
time-dependent network calculation below that critical temperature because the material is 
not completely burned yet. 

In Table 2 we provide the freeze-out temperature and final composition for each of the 
five most abundant nuclei as a function of the initial thermodynamical state and chemical 
composition. As one can see, while the critical temperature for achieving complete NSE was 
the same for all nuclei this was certainly not true for the freeze-out temperature. When the 
initial fuel was composed of carbon and oxygen with no neutron excess, r\ = (rj = 1 — 2Y e ), 
the resulting abundances were largely dominated by 56 Ni. Under the conditions of our 
numerical experiment, 56 Ni froze at T 9 = 3.28, while other Fe-peak isotopes were still reacting 
with residual protons and neutrons. Therefore their freezing temperatures were lower. A 
similar calculation (same initial density and temperature) is shown in the second row of 
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Table 2, but the initial composition was 30% 12 C, 30% 16 and 40% 22 Ne, resulting in a 
moderate neutron excess of i] = 0.0364. This time the dominant nucleus after the adiabatic 
cooling was not 56 Ni but instead were 54 Fe and 58 Ni in accordance with well-known results 
(i.e., Clayton 1968). There is, however, a clear difference in the freezing temperatures of 
these dominant isotopes. 

The third row of Table 2 and Figure 4 summarizes the evolution of a system initially 
composed of pure helium at p = 4 x 10 6 g.cm~ 3 . In this case rj = and the combined 
value of temperature and density in NSE resulted in a larger entropy than in the two earlier 
cases. Now, the period of relaxation to NSE was larger than in the two precedent tests; 
however there was no sign of instability during the integration. Once the NSE was achieved, 
at T = 4.3 x 10 9 K, the composition was again dominated by Fe-peak elements with a lower 
concentration of light particles a and p, (0.66% and 0.74% by mass respectively) than in 
the two precedent cases. For a given initial density the outcome of the adiabatic cooling 
depends on the adopted expansion rate. For very fast expansions the alpha particles have 
some difficulty to react, owing the low-density environment, and the final distribution is rich 
in a particles and intermediate- mass elements (Woosley, Arnett & Clayton 1973). In our 
calculation the expansion rate was not rapid enough to allow the synthesis of intermediate- 
mass elements in appreciable amounts. According to Table 2 the final composition was 
completely dominated by 56 Ni (98%) and Fe-peak elements plus a little amount of unburnt 
helium. 



3.1.1. Influence of screening factors in the equilibrium distribution of nuclei 

As mentioned above, the enhancement factors to the nuclear reaction rates given by 
equation (11) were taken into account in all the realized tests. Generally speaking, their 
inclusion in explosive processes, although important, is not very crucial; its main effect 
being to shift the temporal evolution and the maximum value achieved by temperature. 
Therefore, there is expected to be some changes in the mass fractions of the species in NSE 
when screening factors are incorporated into the scheme, especially at high density. In Figure 
5 there is represented the evolution of temperature during the isochoric combustion of a 50%- 
50% mix of carbon and oxygen at constant density of 4 x 10 9 g.cm^ 3 . As we can see the 
profiles with and without enhancement factors are not very different, although the induction 
time is shorter and the equilibrium temperature is no longer the same: T 9 = 9.51 (no 
screening included) and T 9 = 9.77 (screening included). Again, there is a stable numerical 
behaviour of the system when the equilibrium is approached. 

To analyze the nucleosynthesis during the NSE stage it is more convenient to compare 
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at the same temperature. In doing that not only the comparison among abundances becomes 
more reliable, but a direct check with the results obtained by using the nuclear Saha equation 
with the mass excesses corrected from the coulombic part of the chemical potentials deduced 
from equation (10) is also possible. The results of such comparison for the 14-nuclei network 
are depicted in Figure 6. From that figure it becomes clear that the main effect of screening 
corrections is to shift the distribution of the species towards higher atomic masses and that, 
the larger the density the bigger the change. In fact, the change in the NSE distribution of 
the nuclei depicted in Figure 6 looks rather spectacular. This is in part due to the limited size 
of the 14-nuclei chain; when the same calculation was repeated with the 86-nuclei network 
such differences were spread out and smoothed. For example, at T 9 = 9.32 ,p 9 = 4 the 
mean molecular weight calculated from the a-network was rii = 11.04 g.mol -1 (no screening 
factors) and /ij = 14.35 g.mol -1 (with screening factors) whereas the corresponding values 
for the 86-nuclei network were /Xj = 10.54 g.mol -1 and ^ = 11.87 g.mol -1 respectively. 
The comparison between the resulting mass fractions calculated through the network and 
the nuclear Saha equation also confirms that the adequate way to handle the coulombic 
interactions in the NSE regime is by substracting the electrostatic chemical potentials from 
the nuclear mass excesses (Bravo&Garcia-Senz 1999). 

3.1.2. Iterative refinement and accuracy 

The integration scheme developed in § 2 can be considered as the first iteration step 
of a more general Newton-Raphson based method for solving the full non-linear equations 
governing the evolution of the nuclear system. When using these methods there is often 
an enhancement of the computational efficiency despite that the time invested per model 
is larger. In addition they also have the advantage that a check on accuracy is possible, 
allowing for a more consistent choice of the time-steps (Timmes 1999). Fortunately it is 
straightforward to modify our numerical scheme, with minimum changes, in order to allow 
the iterative refinement of the solution, as described below. The resulting algorithm is 
equivalent to the standard implicit multidimensional Newton-Raphson calculation. For a 
given n + 1 time-step the iterative process leads to a sequence of 'solutions' {Y l n+1 ) k and 
(T n+1 ) h for the i-specie and temperature which, after m iterations, should converge to 
the true values for Y t n+1 and T n+1 . We define a set of functions Gf, which represents the 
discretized form of the chemical equations (1): 

G °[(Y" +1 ) fe , (T n+1 ) k \ = F l [(Y n+1 ) k , (T n+1 f \ At - [(Y l n+1 ) k - Y t n ] k = l...m (13) 
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where (Y n+1 ) fe = {(Y{ l+1 ) k ....{Y™ +l ) k } and F t are the terms on the right of equation (1). A 
similar function representing the energy equation (8) is also defined: 

The Newton- Raphson iterative method demands that (SG^' 1 = — G^ 1 , which allows to 
solve for the corrections AY i: AT to the approximate values (Y t n+1 ) k and (T n+1 ) k found in 
the latest iteration. In this case the Jacobian matrix associated to 8G?' 1 is exactly the same 
as the matrix of the linear system associated to equations (5) and (8). For the first iteration 
we take (F i n+1 ) 1 = Y? and (T n+1 ) 1 = T n so that when m=l the extended Arnett&Truran 
method given by equations (5) and (8) is recovered. For m > 1 an iterative sequence results 
which usually ends when the /c-corrections AY k to (Y™ + ) fe_1 and AT k to (j ln+1 ) k ~ 1 becomes 
negligible or, better, when the G®' 1 functions defined by equations (13) and (14), conveniently 
normalized, go to zero (for example we took |G°|/v/(FiAt) 2 + Alf < 5 x 1(T 6 and a similar 
expression for G 1 ) . Usually convergence is achieved in 3-4 iterations. The criterion to 
choose the time-step was not very restrictive, we allowed for a 1.5% relative variation in 
temperature and 10% in molar fractions. In Figure 7 there is shown the evolution of the 
errors in mass fractions of several species during the isochoric carbon combustion, estimated 
by comparing the size of the corrections given by the first iteration with the algebraic sum of 
the lowest-order corrections given by the remaining iterations. As we can see in that figure 
the first iteration always led to a much larger correction, usually of two orders of magnitude 
or more, than that of the sum of the remaining iterations. Thus, the evolution of the 
mass fractions and temperature can reasonably be described taking only one iteration and 
neglecting second and higher-order corrections. For this particular test, a further relaxation 
of the criterion that restrict the time step led to several negative abundances. Nevertheless, 
the iterative Newton-Raphson scheme would also be in trouble in this case, due to the lack 
of convergence. Therefore, taking only one iteration along as a moderate time-step would 
generally be enough, especially if we are mainly interested in the thermal evolution of the 
system and in the gross features of the nucleosynthesis. For larger time-steps or for detailed 
nucleosynthetic studies it is better to rely on the iterative scheme which has the additional 
advantage that a check of the accuracy is feasible (Timmes 1999). 



-14- 



3.2. Hydrostatic and explosive silicon burning 

Hydrostatic combustion of silicon is a quasi-equilibrium process from the beginning due 
to the important role played by photodisintegration reactions. One goal of this calculation 
is to check the ability of our integration scheme to describe silicon photodissociation at 
constant temperature and density. We also studied explosive silicon burning within an 
expanding medium, until the reactions were totally quenched. In both cases we compared 
our results with those of Timmes et al. (2000), which were calculated under the same 
environmental conditions. Thus, a system composed of 100% 28 Si at T=6 x 10 9 K and 
p = 10 7 g.cm -3 was considered. In the hydrostatic case, the temperature was kept constant, 
therefore the feedback between nuclear kinetics and thermal evolution no longer took place. 
Nevertheless, our scheme was also able to keep track of this limiting case without introducing 
any modification. Instead to rewrite our code to not include AT terms in the matrix we 
have artificially raised the value of the specific heat during the silicon dissociation to keep the 
temperature constant. In the adiabatic expansion test the specific heat was again restored 
to its real value. 

The results of these calculations for the small network are shown in Figures 8,9 and 
10. A comparison between the evolution of the nuclear energy generation rate, provided in 
Figure 8, with respect to that shown in Figure la of Timmes et al. (2000) is satisfactory. 
The main difference is the tiny bump shown around t~ 2 x 10 -4 s, which was probably due to 
the different libraries used to compute the nuclear rates. A comparison of their abundances 
(their Fig. lc) with ours (Fig. 9) also shows a good agreement. 

As a second check a piece of pure silicon was enforced to follow an adiabatic expansion 
in the hydrodynamic timescale as shown in equation (12). In this case, the temperature 
evolution was calculated self-consistently from the network. The results of this calculation 
for our a— chain are summarized in Figures 8 and 10. These figures can be compared 
with Figures 2a,2c in Timmes et al. (2000), in which the temperature is evolved with 
an approximate analytical formula. The evolution of the nuclear energy rate follows a close 
path in both calculations. The comparison between abundances is also quite satisfactory, 
although there are divergences above t=0.1 s, which is about the hydrodynamical time thd- 
These differences probably arise from the slightly different path followed by temperature 
in both calculations. We note that Timmes et al. (2000) only included the contribution 
of the radiation to the specific heat thus, for a given density, we found a higher value of 
temperature. In particular, our final silicon group abundance (up to titanium) was about 
five times higher than theirs. On the whole its composition is characteristic of an a— rich 
freeze-out with a large abundance of 56 Ni (92% in mass) followed by 4 He (7%). 
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3.3. An example of isobaric combustion: The propagation of a subsonic 

nuclear flame 

Thermonuclear flame propagation is a subsonic mode of energy propagation which takes 
place under nearly isobaric conditions. Flames could appear due to the combination of 
sharp thermal gradients, high conductivity and potentially explosive fuels such as helium 
or carbon and oxygen. These requirements may be fulfilled in the interior of those massive 
white dwarfs that are reactivated by the accretion from a companion star, giving rise to a 
Type la supernova. This calculation is more complex than those carried out in the previous 
subsections because both, the temporal evolution and the spatial structure of the flame have 
to be solved. Therefore our nuclear subroutine was coupled to the diffusion equation within 
an initially uniform grid of shells representing a physical system with planar geometry and 
constant density. Given the adequate initial conditions (see below), the advance of the flame 
was calculated directly from the nuclear network of 86 isotopes along with the isobaricity 
assumption (Timmes & Woosley 1992). A diffusive term was added to the energy equation 
so that the heat is transported from the burned to the unburned zone, mainly by electronic 
conduction. Thus, we take, 

dQ = S dif dt=-V ■ (ffVT) dt (15) 

in equation (9) and set P ext = constant in the same equation in order to keep the pressure 
constant during the flame propagation. The conductivity, a, was that of Khokhlov, Oran 
and Wheeler (1997) which includes electronic and radiative contributions to the opacity. As 
the dependence of conductivity on temperature is much lower than it is in nuclear reactions 
this new term, equation (15), can be added to the independent terms in the corresponding 
discretized version of equation (9) . In order to build a nuclear flame the left side of a 
tube containing equal parts of carbon and oxygen was incinerated at constant pressure until 
NSE was achieved. After a brief transitory period the heat exchange between the burned 
the and cold material gave rise to a steady planar flame that moved at a constant velocity 
through the fuel with constant density Pf ue i = 1-26 x 10 8 g.cm -3 . Behind the flame front the 
original fuel was rapidly transformed into intermediate-mass elements, especially silicon; the 
temperature rose to 6.1 x 10 9 K; and the density dropped to p = 7.6 x 10 7 g.cm -3 to keep the 
pressure constant. Under these conditions the ashes were in QSE and began to slowly relax 
to complete equilibrium, as can be seen in Figure 11. During the total elapsed time (about 
1CT 5 s) the nuclear system did not show any sign of instability. The stationary velocity of 
the flame front was v// ame = 4.8 km.s -1 , which is about 50% higher than that obtained by 
Timmes & Woosley (1992) for the same density. Such discrepancy is well within the factor 
two error usually attributed to the input physics (conductivities, nuclear libraries and EOS). 
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4. Conclusions 

Current standard integration methods for chemical equations encounter some difficulty 
when handling high temperature combustion near equilibrium. The reason for this is the 
extreme sensitivity of nuclear rates to the temperature, which makes it difficult to reach a 
time-independent situation (equilibrium) when using a time-dependent solver for equations 
that describe the evolution of molar fractions. Miiller (1986) pointed out that the solution 
is to implicitly couple the equations that give abundance changes with those that describe 
energy or entropy evolution. A practical drawback of such a procedure is that the derivatives 
of the nuclear reaction rates have to be worked out at each iteration. Nowadays, however, 
very complete libraries are available, which detail nuclear reaction rates for many species 
with simple formulae (i.e. Rauscher & Thielemann 2000). From these, then, derivatives can 
be evaluated without causing excessive computational overload. 

The main purpose of this study was to build and check an integration method, inspired 
by Miiller idea, that would be flexible enough to accommodate a variety of combustion 
modes involving quasi/total nuclear equilibrium at some point along the path of its evolu- 
tion. Screening enhancement factors were also incorporated to our numerical scheme. Good 
astrophysical examples of the potential applicability of this algorithm would be for nuclear 
flames and detonation waves in Type la supernovae and for silicon burning in Type II su- 
pernovae respectively. 

In spite of the simplicity of the proposed numerical algorithm, described in §2, we 
found that it was able to pass many tests related to high-temperature combustion, §3. In 
each of them stability was always preserved, even in the NSE stage, without restricting 
the time-step. The results also matched well those obtained in similar conditions by other 
authors using different integration schemes. A comparison of the numerical efficiency of the 
small chain (of 14 nuclei) and the moderate-sized one (of 86 nuclei) showed that the time 
overload introduced by thermal coupling was under 15%, and 2% for the small and large 
networks respectively. These factors arise mainly from the calculation of the derivatives of 
the nuclear rates rather than from the inversion of the matrix. Recently, the 14-isotope 
a-chain with thermal coupling was incorporated to a multidimensional SPH code designed 
to simulate Type la supernova explosions. In that situation, the average time degradation 
is far below the aforementioned 15% as the time expended per model is largely controlled 
by hydrodynamics, especially when one is searching for neighbours or calculating gravity. 

The FORTRAN subroutines that implement the 14 a-chain and the 86-isotope network 
presented in this paper are made available by the authors to interested readers upon request. 
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Fig. 1. — Sketch of the nuclear network used in the tests. It consists of 86 nuclei (including 
p, n and a particles), linked trough (p,7) -thin solid arrows-, (11,7) -small dotted arrows-, 
(a/y) -thick solid arrows-, (a,p) -long dashed arrows-, and (a,n) -long dotted arrows-, as 
suggested by the plot. The network is organized into groups around an a-chain of 14 nuclei 
from 4 He to 60 Zn. One of these groups, between 20 Ne and 24 Mg, is shown in the figure. 
Triple alpha reaction and binary reactions of carbon and oxygen were also included. 



-20- 




Model 



Fig. 2. — Evolution at constant density, p = 10 9 g.cm~ 3 of a 50% mix of carbon and oxygen 
calculated without implicit thermal coupling, the initial temperature was T = 10 9 K. As 
soon as the combustion enters in the QSE zone and NSE plateau it becomes unstable (86- 
isotope network). 



200 400 600 800 1000 1200 1400 

Model 

Fig. 3. — The same as in Figure 2, this time with the thermal coupling turned on. The 
system remains stable all the time. An adiabatic expansion was imposed above model 600 
until all the nuclear reactions quenched (86-isotope network). 
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Fig. 4. — Isochoric combustion of pure helium (at p = 4 x 10 6 g.cm" 3 ) followed by an adia- 
batic expansion in the hydrodynamic time. Owing to the lower ignition density it took more 
time for the combustion to reach complete equilibrium at T = 4.23 x 10 9 K. Nevertheless, 
the implicit thermal coupling again avoids any instability to grow. The resulting nucleosyn- 
thesis after the expansion is very rich in 56 Ni, with appreciable amounts of helium and minor 
concentration of a-elements such as 40 Ca (86-isotope network). 
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Fig. 5. — Influence of the screening corrections on the evolution of temperature during the 
isochoric combustion of a 50% mix of carbon and oxygen at p = 4 x 10 9 g.cm -3 . Inclusion of 
screening factors (dashed line) led to a higher temperature value during the NSE. (14-isotope 
chain). 
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Fig. 6. — Abundance of 4 He, 28 Si and 56 Ni at several densities once the NSE has been 
achieved. For a given density and temperature the isotopic composition at NSE sensitively 
depends on the inclusion of the screening corrections to the nuclear reaction rates. Contin- 
uum lines depicts the abundances when no screening factors are taken into account. Dashed 
lines are for the screening factors given by equation (11). A comparison with the mass frac- 
tions obtained by solving the nuclear Saha equation for the same number of species (14) is 
also provided (triangles and diamonds). In this case the electrostatic interactions were taken 
into account by substracting the chemical potential of the nuclei from their nuclear mass 
excess. There is a good agreement between both calculations. As expected, the effect of the 
electrostatic interactions increases with density. (14-isotope chain). 
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Fig. 7.— Evolution of AX//X; (solid lines), and J2 k AX i/ X i> k = 2,m (dashed lines) 
corresponding to He 4 , C 12 and Si 28 and the equivalent expressions for temperature, during 
the isochoric combustion of a 50% mix of carbon and oxygen at p = 10 9 g.cm -3 . Horitzontal 
axis is temperature, which goes from T° = 10 9 K to T = 8.6 x 10 9 just before the relaxation 
to NSE. The error during the purely linear calculation (one iteration) is approximately given 
by the sum of the higher-order corrections (dashed lines) calculated through the iterative 
Newton-Raphson scheme described in §3.1.2. The time-step was chosen to allow a relative 
variation of 1.5% and 10% in temperature and mass-fractions of the more abundant species 
respectively (86-isotope chain). 
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Fig. 8. — Evolution of the nuclear energy generation rate (absolute value) during endoener- 
getic silicon combustion: at constant density p = 10 7 g.cm -3 and temperature T = 6 x 10 9 K 
(dashed line), and during the explosive combustion (solid line). For the last case the nuclear 
energy generation rate becomes positive after 1CT 3 s, (14-isotope chain). 
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Fig. 9. — The same as Figure 8 but for abundances in the hydrostatic calculation (14-isotope 
chain) . 
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Fig. 10. — The same as in Figure 9 but for the hydro dynamic silicon combustion. Final 
products are typical of an a-rich freeze out (14-isotope chain). 
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Fig. 11. — Temperature and chemical profiles corresponding to a nuclear flame moving to the 
right through a medium of density p = 1.26 x 10 8 g.cm~ 3 composed of equal parts carbon and 
oxygen. Again, the feedback between nuclear rates and temperature was adequately handled 
during the integration, and the relaxation to NSE was smooth (86-isotope network). 
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Table 1. Abundances during the adiabatic expansion from NSE at 
po = 10 9 g.cuT 3 starting from the NSE temperature T = 8.03 10 9 K. 



Isot. 


T 9 : 8.03 


4 


3.0 


2.0 


0.5 


Isot. 


T 9 : 8.03 


4.0 


3.0 


2.0 


0.5 


n 


5.7 


-5) 


4.8(-12) 


3.0(-16) 


8.6(-25) 
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38 Ca 
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-7) 


3.1(-14) 
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l.l(-8) 
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P 
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1.7(-3) 
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39 Ca 


2.3 
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-1) 
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4.4(-3) 


40 Ca 


5.6 


-3) 


6.6(-8) 


1.3(-7) 
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12 C 


2.5 


-5) 
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4.1 (-7) 


5.1(-7) 


41 Ca 


1.3 


-3) 


1.7(-11) 


1.9(-12) 


1.5(-18) 


1.3(-21) 


16 


5.3 
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1.7(-9) 


3.2(-9) 


1.1 (-9) 


3.7(-10) 


42 Ca 
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-4) 


6.4(-13) 


7.6(-15) 
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6.8(-27) 


20 Ne 


3 
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3.3(-ll) 
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1.6(-11) 


4.7(-12) 


41 Sc 
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2(-44) 
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6.6(-13) 
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3(-19) 
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5 
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21 Na 


4 
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42 Ti 
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8(-10) 
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9.4(-15) 
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2 (-40) 


43 Ti 
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2.6(-13) 
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8.3(-12) 


8.3(-12) 
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44 Ti 
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6(-8) 


7.9(-8) 
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5.9(-8) 


45 Ti 
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46 Ti 
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24 Mg 
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46 Cr 
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2(-13) 
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26 A1 
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4 
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6 
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1.4(-15) 
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2.8(-17) 
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9.63(-l) 
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8 
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5.5 
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4.5 
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58 Ni 
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-1) 
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7.9(-10) 
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-5) 
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37 Ar 
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-3) 


1(-12) 
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5 


-4) 


6.2(-5) 


1.4(-4) 


6.2(-3) 


6.8(-3) 


38 Ar 


1.1 


-3) 


1.9(-13) 


5.4(-15) 


8(-22) 


6.5(-30) 


59 Cu 


5.4 


-3) 


6(-4) 


2.4(-4) 


5.8(-7) 


1.7(-7) 


37 K 


2.8 


-5) 


6.8(-13) 


3(-12) 


1.2(-11) 


7.1(-11) 


58 Zn 


1.2 


-8) 


6(-10) 


4.5(-10) 


1.4(-8) 


2.7(-9) 


38 K 


2.3 


-4) 


2.3(-12) 


1.6(-13) 


7(-18) 


l.l(-22) 


59 Zn 


1.8 


-6) 


8.7(-8) 


6.7(-7) 


2.8(-3) 


1.02(-2) 


39 K 


2.3 


-4) 


3.1(-11) 


7.9(-ll) 


8.5(-13) 


1.1(-13) 


60 Zn 


1.3 


-4) 


1.5(-4) 


1.7(-3) 


4.9(-3) 


4.9(-3) 



-31 - 



-32 - 



Table 2. Five more abundant nuclei at NSE and at the freeze-out temperature for the 
models considered in Sec. 3.1. The NSE temperature and mass fractions of these nuclei are 
given in the left side of the Table, as a function of the initial composition of the fuel. The 
approximate freeze out temperature and the corresponding abundances are also given in 
the right part of the Table. A specie is considered burned once its mass fraction differs 
from its final stable value in less than 2%. Temperatures are in units of 10 9 K. 



Fuel 






NSE 










FREEZE 


OUT 








Ash 


54 Fe 


4 He 


58 Ni 


55 Co 


57 Ni 


56 Ni 


59 Zn 


57 Ni 


58 Cu 


60 Zn 


CO 


Xu 


2.7(-l) 


1.95(-1) 


1.14(-1) 


9.55(-2) 


5.46(-2) 


9.63(-l) 


1.02(-2) 


1.01(-2) 


6.82(-3) 


4.93(-3) 




T 9 


8.03 


8.03 


8.03 


8.03 


8.03 


3.20 


1.31 


1.37 


1.35 


2.10 




Ash 


54 Fe 


4 He 


50 Cr 


58 N ; 


34 g 


54 Fe 


58 Ni 


50 Cr 


4 He 


46 Ti 


CONc 


A"fc 


4.1(-1) 


1.99(-1) 


1.59(-1) 


1.04(-1) 


3.08(-2) 


7.25(-l) 


2.12(-1) 


5.58(-2) 


9.49(-3) 


3.11(-3) 




T 9 


8.51 


8.51 


8.51 


8.51 


8.51 


5.52 


4.77 


4.75 


4.45 


4.63 




Ash 


56 Ni 


85 Co 


54 Fe 


w Ni 


58 Ni 


56 Ni 


57 Ni 


59 Zn 


58 Cu 


4 He 


He 


x k 


6.2(-l) 


l.l(-l) 


9.3(-2) 


5.0(-2) 


3.8(-2) 


9.7(-l) 


8.2(-3) 


6.5(-3) 


4.3(-3) 


1.7(-3) 




T 9 


4.30 


4.32 


4.32 


4.32 


4.32 


2.95 


1.14 


1.11 


1.23 


2.43 



